BayesExamples
Kole Norberg
4/19/2021
R Markdown
First up is to load in the necessary packages. You may need to install them. brms uses stan which can be a pain to load. It might ask you if you want to compile. If I’m having trouble, sometimes I just say no and it works. Also, start with a clean R environment when you install and load up brms or anything else with stan in it. There are some libraries that can get in the way. I don’t really know the source of the difficulty. I just know sometimes it installs quickly and easily and sometimes it doesn’t. You may need to google around a bit if you’re having trouble.
#install.packages("brms")
#install.packages("BayesFactor")
#install.packages("simstandard")
#install.packages("tidyverse")
library(brms)## Loading required package: Rcpp## Loading 'brms' package (version 2.15.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').##
## Attaching package: 'brms'## The following object is masked from 'package:stats':
##
## arlibrary(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.3.2 ✓ purrr 0.3.4
## ✓ tibble 3.0.4 ✓ dplyr 1.0.2
## ✓ tidyr 1.1.2 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.0## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(simstandard)
library(BayesFactor)## Loading required package: coda## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack## ************
## Welcome to BayesFactor 0.9.12-4.2. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
##
## Type BFManual() to open the manual.
## ************##Data simulation I’ve created four data sets. Two have 40 observations and two have 400. Priors will affect the results differently depending on the number of observations. Each set (small and large) also have a case in which the data will come out to be sig (p<.05) and a case where they will not (p>.05). Note that in the sig case, the data is standardized (the intercept will therefore be close to 0). I’ve posted the histograms so you can see how the distributions look and have computed the correlations between the x (age) and y (vocabulary knowledge) variables to verify that the relationship in the significant set is greater than in the non-sig. set. Using the same random see (set.seed()) will usually ensure you get the same results as I do.
set.seed(4321)
df.ns.small <- tibble::tibble(vk = rnorm(40,200,50), age=rnorm(40,18,4))
summary(df.ns.small)## vk age
## Min. :120.6 Min. :12.77
## 1st Qu.:178.1 1st Qu.:17.48
## Median :201.1 Median :18.74
## Mean :203.0 Mean :18.95
## 3rd Qu.:233.8 3rd Qu.:21.12
## Max. :280.5 Max. :26.32hist(df.ns.small$age)
hist(df.ns.small$vk)
cor(df.ns.small$age,df.ns.small$vk) #.02## [1] 0.01946101df.s.small <- simstandard::sim_standardized("age ~~ 0.75 * vk", n = 40)
summary(df.s.small)## age vk
## Min. :-2.00079 Min. :-1.52273
## 1st Qu.:-0.81692 1st Qu.:-0.50523
## Median : 0.01170 Median : 0.01878
## Mean :-0.06033 Mean : 0.08363
## 3rd Qu.: 0.82217 3rd Qu.: 0.82717
## Max. : 2.23053 Max. : 1.95663hist(df.s.small$age)
hist(df.s.small$vk)
cor(df.s.small$age,df.s.small$vk) #.64## [1] 0.6385461set.seed(4321)
df.ns.large <- tibble::tibble(vk = rnorm(400,200,50), age=rnorm(400,18,4))
summary(df.ns.large)## vk age
## Min. : 77.5 Min. : 5.898
## 1st Qu.:170.5 1st Qu.:15.148
## Median :202.4 Median :17.909
## Mean :203.0 Mean :17.804
## 3rd Qu.:236.2 3rd Qu.:20.478
## Max. :327.0 Max. :34.077hist(df.ns.large$age)
hist(df.ns.large$vk)
cor(df.ns.large$age,df.ns.large$vk)#.02## [1] 0.02486447df.s.large <- simstandard::sim_standardized("age ~~ 0.75 * vk", n = 400)
summary(df.s.large)## age vk
## Min. :-2.81939 Min. :-2.84270
## 1st Qu.:-0.56397 1st Qu.:-0.63119
## Median : 0.02679 Median :-0.02617
## Mean : 0.01714 Mean : 0.01719
## 3rd Qu.: 0.66817 3rd Qu.: 0.60855
## Max. : 2.92719 Max. : 3.52546hist(df.s.large$age)
hist(df.s.large$vk)
cor(df.s.large$age,df.s.large$vk)#.71## [1] 0.7136371##Setting priors I’m going to do everything in brms because I think it’s the easiest to work with all things considered. stan_lm in the rstanarm library is also good.
I’m setting the priors below to be either diffuse (which here is still somewhat informative), weak, or informative based on a normal distribution. If you’re doing logistic regression, you would actually probably want to set the priors on a beta distribution (not the same beta as the beta coefficient) which can be a bit tricky, so I suggest watching a few videos on the beta distribution if you plan to do logistic regression with brms. The cauchy distribution is also commonly used with continuous data. I’m using the normal here because it’s more familiar, but you’ll want to look up the distribution you use to make sure it reflects your prior beliefs.
For the normal, I set the mean and sd. The mean is always the same in this example; it’s the sd that determines how informative the prior is. Class refers to either the variable (b - think beta coefficient), the intercept, sigma (error), or random interecepts and slopes if you’re using mixed models. Coef is the name of the variable (case sensitive). The priors I’m setting up here suggest a belief that for every increase in age by one year, vocabulary knowledge increased by 5 but that the mean (intercept) for vocabulary knowledge was 200 in the case of raw data or 0 in the case of standardized data.
Getting a bayes factor requires setting reasonable priors. The way you set the prior can heavily bias your bayes factor. The best way to handle this is to do a sensitivity analysis to see how your bayes factor reacts to different levels of priors. Here we will also look at how it reacts to priors that suggest there is no relationship (null).
Although we’re primarily interested in the beta coefficient on the variable, setting a prior for the intercept can also be useful. Here I’m using the student_t because it’s normally used for the intercept in default cases. The first number you get is the df, then the mean, then the sd. Because we’re mostly interested in the deviation from the intercept, I think it’s ok to go with the default prior from brms and have only set the prior on the intercept once as an example.
#Priors for an effect, not standardized
diffuse <- c(set_prior("normal(5,100)", class = "b", coef = "age"),
set_prior("student_t(3,200,2.5)", class="Intercept"))
weak <- set_prior("normal(5,10)", class = "b", coef = "age")
informative <- set_prior("normal(5,1)", class = "b", coef = "age")
hyperinformative <- set_prior("normal(5,.1)", class = "b", coef = "age")
#Priors for an effect, standardized. I only need to do thie for the diffuse prior because it has a prior on the intercept which will change when standardized
diffuse.s <- c(set_prior("normal(5,100)", class = "b", coef = "age"),
set_prior("student_t(3,0,2.5)", class="Intercept"))
#Null Priors
diffuseN <- c(set_prior("normal(0,100)", class = "b", coef = "age"),
set_prior("student_t(3,200,2.5)", class="Intercept"))
diffuse.sN <- c(set_prior("normal(0,100)", class = "b", coef = "age"),
set_prior("student_t(3,0,2.5)", class="Intercept"))
weakN <- set_prior("normal(0,10)", class = "b", coef = "age")
informativeN <- set_prior("normal(0,1)", class = "b", coef = "age")
hyperinformativeN <- set_prior("normal(0,.1)", class = "b", coef = "age")##Running brms brms is set up just like the lm with a few other parameters. family=gaussian() says that this is a continuous variable and we are using a normal distribution. If this were logistic, I would use family=bernoulli (not binomial). The math for these distributions is different but for inferential purposes, they produce the same results. Brms prefers bernoulli.
Chains refers to the markov chains. 4 is plenty. Explaining this is beyond the scope of this document. Basically, the model starts in 4 places on the distribution and checks to see if it ends up in the same spot each time. rhat=1 in the summary output if this works.
iter refers to the number of iterations in each chain. 2000 is usually enough, but if you find that in the summary output, rhat does not = 1, you should increase your iterations as you don’t have a good fit. rhat = 1.01 is not good enough; you want 1. 2000 is the default, so I have it set once in ns.s.flat so you can see how to do it, but I don’t set it again.
cores allows you to use your computers processing power to speed up the computations. You can use as many cores as there are chains provided your computer has that many. Each core will run a chain. I have two cores on my computer, so if I set cores to 2, the entire process takes ~half the time.
sample_prior = TRUE allows us to plot our priors later. It’s just telling brms not to trash this data. If you’re having trouble understanding your prior, plotting it is the easiest way to get a sense of what is happening. This only works though if you set priors; it doesn’t work with brms default priors.
save_all_pars = save_pars(all=TRUE) keeps the appropriate values on hand for calculating the bayes factor.
I’m running six models below. the first is the basic linear model, then a brms model with flat priors. Flat priors are set by brms. They are uninformative and improper (don’t sum to 1) priors and consist basically of placing equal probability on all real numbers. brms calls this “flat priors over the reals.” We then move on to the diffuse, weakly informative, informative, and hyperinformative models. The next output consists of the prior summaries followed by the results summaries. Results summaries are similar to lm summaries with a few small differences that we’ll go over in class. The final bit are the plots of the brms models. Note that the flat model has no prior distribution plotted.
#Basic linear model
mod.lm.ns.small <- lm(vk~age,data=df.ns.small)
#brms with flat priors (brms sets a prior here that is flat but over the reals meaning only covers the real numbers)
ns.s.flat <- brm(vk~age,data=df.ns.small,family=gaussian(),chains=4,iter=2000,cores=2,sample_prior = TRUE,save_pars=save_pars(all=TRUE))## Compiling Stan program...## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1## Start sampling#brms with minimal priors
ns.s.diffuse <- brm(vk~age,data=df.ns.small,family=gaussian(),chains=4,cores=2,prior=diffuse,sample_prior = TRUE,save_pars=save_pars(all=TRUE))## Compiling Stan program...
## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1## Start sampling#brms with weak priors
ns.s.weak <- brm(vk~age,data=df.ns.small,family=gaussian(),chains=4,cores=2,prior=weak,sample_prior = TRUE,save_pars=save_pars(all=TRUE))## Compiling Stan program...
## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1## Start sampling#brms with informative priors
ns.s.inf <- brm(vk~age,data=df.ns.small,family=gaussian(),chains=4,cores=2,prior=informative,sample_prior = TRUE,save_pars=save_pars(all=TRUE))## Compiling Stan program...
## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1## Start sampling#brms with hyper informative priors
ns.s.hyper <- brm(vk~age,data=df.ns.small,family=gaussian(),chains=4,cores=2,prior=hyperinformative,sample_prior = TRUE,save_pars=save_pars(all=TRUE))## Compiling Stan program...
## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1## Start samplingprior_summary(ns.s.flat)## prior class coef group resp dpar nlpar bound
## (flat) b
## (flat) b age
## student_t(3, 201.1, 44.9) Intercept
## student_t(3, 0, 44.9) sigma
## source
## default
## (vectorized)
## default
## defaultprior_summary(ns.s.diffuse)## prior class coef group resp dpar nlpar bound source
## (flat) b default
## normal(5,100) b age user
## student_t(3,200,2.5) Intercept user
## student_t(3, 0, 44.9) sigma defaultprior_summary(ns.s.weak)## prior class coef group resp dpar nlpar bound source
## (flat) b default
## normal(5,10) b age user
## student_t(3, 201.1, 44.9) Intercept default
## student_t(3, 0, 44.9) sigma defaultprior_summary(ns.s.inf)## prior class coef group resp dpar nlpar bound source
## (flat) b default
## normal(5,1) b age user
## student_t(3, 201.1, 44.9) Intercept default
## student_t(3, 0, 44.9) sigma defaultprior_summary(ns.s.hyper)## prior class coef group resp dpar nlpar bound source
## (flat) b default
## normal(5,.1) b age user
## student_t(3, 201.1, 44.9) Intercept default
## student_t(3, 0, 44.9) sigma defaultsummary(mod.lm.ns.small)##
## Call:
## lm(formula = vk ~ age, data = df.ns.small)
##
## Residuals:
## Min 1Q Median 3Q Max
## -83.225 -24.760 -1.386 31.988 77.630
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 198.0078 41.9681 4.718 3.19e-05 ***
## age 0.2625 2.1876 0.120 0.905
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 41.66 on 38 degrees of freedom
## Multiple R-squared: 0.0003787, Adjusted R-squared: -0.02593
## F-statistic: 0.0144 on 1 and 38 DF, p-value: 0.9051summary(ns.s.flat)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vk ~ age
## Data: df.ns.small (Number of observations: 40)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 197.71 43.03 112.18 278.62 1.00 3990 2869
## age 0.28 2.23 -4.05 4.78 1.00 3977 2764
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 42.42 4.84 34.27 53.20 1.00 3365 2795
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).summary(ns.s.diffuse)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vk ~ age
## Data: df.ns.small (Number of observations: 40)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 195.90 41.70 115.36 280.19 1.00 3572 2728
## age 0.24 2.19 -4.16 4.49 1.00 3540 2648
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 42.06 4.75 33.76 52.64 1.00 3884 2413
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).summary(ns.s.weak)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vk ~ age
## Data: df.ns.small (Number of observations: 40)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 193.37 41.23 111.02 274.01 1.00 3468 2873
## age 0.51 2.15 -3.70 4.80 1.00 3463 2891
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 42.47 5.00 34.15 53.52 1.00 3253 2802
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).summary(ns.s.inf)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vk ~ age
## Data: df.ns.small (Number of observations: 40)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 123.04 18.65 86.10 158.42 1.00 3653 2697
## age 4.22 0.92 2.43 6.04 1.00 3649 3215
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 43.77 4.96 35.19 54.51 1.00 3347 2708
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).summary(ns.s.hyper)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vk ~ age
## Data: df.ns.small (Number of observations: 40)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 108.40 7.04 94.64 122.18 1.00 3776 3111
## age 4.99 0.10 4.80 5.19 1.00 3664 2631
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 44.36 5.28 35.55 55.95 1.00 3664 2753
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(hypothesis(ns.s.flat, "age = 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Flat Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)
plot(hypothesis(ns.s.diffuse, "age = 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Diffuse Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)
plot(hypothesis(ns.s.weak, "age = 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Weak Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)
plot(hypothesis(ns.s.inf, "age = 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Informative Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)
plot(hypothesis(ns.s.hyper, "age = 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Hyper Informative Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)
The below graphs just put the diffuse, weak, and informative models together into one graph so you can see how the prior distribution affects the posterior.
posterior1 <- posterior_samples(ns.s.diffuse, pars = "b_age")[,c(1,2)]
posterior2 <- posterior_samples(ns.s.weak, pars = "b_age")[,c(1,2)]
posterior3 <- posterior_samples(ns.s.inf, pars = "b_age")[,c(1,2)]
posterior1.2.3 <- bind_rows("prior 1" = pivot_longer(posterior1,c(prior_b_age,b_age)),
"prior 2" = pivot_longer(posterior2,c(prior_b_age,b_age)),
"prior 3" = pivot_longer(posterior3,c(prior_b_age,b_age)),
.id = "id")
modelLME <- lm(vk ~ 1 + age, data = df.ns.small)
ggplot(data = posterior1.2.3,
mapping = aes(x = value,
fill = id,
colour = name,
linetype = name,
alpha = name
)) +
geom_density(size = 1.2)+
geom_vline(xintercept = summary(modelLME)$coefficients["age", "Estimate"], #add the frequentist solution too
size = .8, linetype = 2, col = "black")+
scale_x_continuous(limits = c(-20, 20))+
coord_cartesian(ylim = c(0, .5))+
scale_fill_manual(name = "Densities",
values = c("Yellow","darkred","blue" ),
labels = c("diffuse ~ N(5,100) prior",
"weak ~ N(5,10) prior",
"informative ~ N(5, 1) prior") )+
scale_colour_manual(name = 'Posterior/Prior',
values = c("black","red"),
labels = c("posterior", "prior"))+
scale_linetype_manual(name ='Posterior/Prior',
values = c("solid","dotted"),
labels = c("posterior", "prior"))+
scale_alpha_discrete(name = 'Posterior/Prior',
range = c(.7,.3),
labels = c("posterior", "prior"))+
annotate(geom = "text",
x = 0.45, y = -.13,
label = "LME estimate: 0.804",
col = "black",
family = theme_get()$text[["family"]],
size = theme_get()$text[["size"]]/3.5,
fontface="italic")+
labs(title = expression("Influence of Priors on Posterior"),
subtitle = "3 different densities of priors and posteriors and the LME estimate")+
theme_bw()## Warning: Using alpha for a discrete variable is not advised.## Warning: Removed 3650 rows containing non-finite values (stat_density).
We will do all of the above again, except this time for the data set we know is significant. I’m now using include = FALSE in the header in order to suppress the long progress output within the markdown file.
summary(mod.lm.s.small)##
## Call:
## lm(formula = vk ~ age, data = df.s.small)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.76709 -0.37218 -0.07462 0.45164 1.17837
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1155 0.1060 1.090 0.283
## age 0.5284 0.1033 5.115 9.27e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6691 on 38 degrees of freedom
## Multiple R-squared: 0.4077, Adjusted R-squared: 0.3922
## F-statistic: 26.16 on 1 and 38 DF, p-value: 9.271e-06summary(s.s.flat)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vk ~ age
## Data: df.s.small (Number of observations: 40)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.11 0.11 -0.10 0.33 1.00 3278 2563
## age 0.53 0.11 0.33 0.74 1.00 3162 2944
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.69 0.08 0.55 0.88 1.00 3118 2724
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).summary(s.s.diffuse)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vk ~ age
## Data: df.s.small (Number of observations: 40)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.12 0.11 -0.10 0.33 1.00 3365 2981
## age 0.53 0.10 0.32 0.73 1.00 3523 2835
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.69 0.08 0.55 0.88 1.00 3036 2615
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).summary(s.s.weak)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vk ~ age
## Data: df.s.small (Number of observations: 40)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.11 0.11 -0.11 0.33 1.00 3941 3054
## age 0.53 0.11 0.32 0.75 1.00 3726 2708
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.69 0.08 0.55 0.88 1.00 3615 2723
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).summary(s.s.inf)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vk ~ age
## Data: df.s.small (Number of observations: 40)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.12 0.11 -0.11 0.34 1.00 2868 2509
## age 0.58 0.11 0.37 0.80 1.00 3744 2717
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.70 0.08 0.56 0.88 1.00 3481 3029
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).summary(s.s.hyper)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vk ~ age
## Data: df.s.small (Number of observations: 40)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.38 0.70 -1.00 1.77 1.00 3712 2791
## age 4.91 0.10 4.71 5.11 1.00 3254 2745
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 4.60 0.52 3.69 5.79 1.01 3382 2468
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(hypothesis(s.s.flat, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Flat Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)
plot(hypothesis(s.s.diffuse, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Diffuse Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)
plot(hypothesis(s.s.weak, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Weak Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)
plot(hypothesis(s.s.inf, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Informative Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)
plot(hypothesis(s.s.hyper, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Hyper Informative Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)
posterior1 <- posterior_samples(s.s.diffuse, pars = "b_age")[,c(1,2)]
posterior2 <- posterior_samples(s.s.weak, pars = "b_age")[,c(1,2)]
posterior3 <- posterior_samples(s.s.inf, pars = "b_age")[,c(1,2)]
posterior1.2.3 <- bind_rows("prior 1" = pivot_longer(posterior1,c(prior_b_age,b_age)),
"prior 2" = pivot_longer(posterior2,c(prior_b_age,b_age)),
"prior 3" = pivot_longer(posterior3,c(prior_b_age,b_age)),
.id = "id")
modelLME <- lm(vk ~ 1 + age, data = df.s.small)
ggplot(data = posterior1.2.3,
mapping = aes(x = value,
fill = id,
colour = name,
linetype = name,
alpha = name
)) +
geom_density(size = 1.2)+
geom_vline(xintercept = summary(modelLME)$coefficients["age", "Estimate"], #add the frequentist solution too
size = .8, linetype = 2, col = "black")+
scale_x_continuous(limits = c(-20, 20))+
coord_cartesian(ylim = c(0, .5))+
scale_fill_manual(name = "Densities",
values = c("Yellow","darkred","blue" ),
labels = c("diffuse ~ N(5,100) prior",
"weak ~ N(5,10) prior",
"informative ~ N(5, 1) prior") )+
scale_colour_manual(name = 'Posterior/Prior',
values = c("black","red"),
labels = c("posterior", "prior"))+
scale_linetype_manual(name ='Posterior/Prior',
values = c("solid","dotted"),
labels = c("posterior", "prior"))+
scale_alpha_discrete(name = 'Posterior/Prior',
range = c(.7,.3),
labels = c("posterior", "prior"))+
annotate(geom = "text",
x = 0.45, y = -.13,
label = "LME estimate: 0.804",
col = "black",
family = theme_get()$text[["family"]],
size = theme_get()$text[["size"]]/3.5,
fontface="italic")+
labs(title = expression("Influence of Priors on Posterior"),
subtitle = "3 different densities of priors and posteriors and the LME estimate")+
theme_bw()## Warning: Using alpha for a discrete variable is not advised.## Warning: Removed 3670 rows containing non-finite values (stat_density).
We’re repeating the above again for large not significant data set
summary(mod.lm.ns.large)##
## Call:
## lm(formula = vk ~ age, data = df.ns.large)
##
## Residuals:
## Min 1Q Median 3Q Max
## -125.933 -31.983 -0.136 33.682 122.293
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 197.8222 10.7598 18.385 <2e-16 ***
## age 0.2924 0.5893 0.496 0.62
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 47.66 on 398 degrees of freedom
## Multiple R-squared: 0.0006182, Adjusted R-squared: -0.001893
## F-statistic: 0.2462 on 1 and 398 DF, p-value: 0.62summary(l.ns.flat)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vk ~ age
## Data: df.ns.large (Number of observations: 400)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 197.59 10.69 176.94 218.56 1.00 4103 3106
## age 0.31 0.59 -0.84 1.44 1.00 4172 3159
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 47.77 1.72 44.46 51.17 1.00 4172 2834
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).summary(l.ns.diffuse)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vk ~ age
## Data: df.ns.large (Number of observations: 400)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 196.38 10.81 175.44 216.96 1.00 4021 2807
## age 0.30 0.60 -0.85 1.48 1.00 4000 2827
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 47.76 1.66 44.65 51.22 1.00 3886 2895
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).summary(l.ns.weak)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vk ~ age
## Data: df.ns.large (Number of observations: 400)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 197.80 10.79 176.18 218.85 1.00 4463 3018
## age 0.29 0.59 -0.88 1.46 1.00 4462 3154
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 47.76 1.71 44.55 51.26 1.00 4132 2803
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).summary(l.ns.inf)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vk ~ age
## Data: df.ns.large (Number of observations: 400)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 176.02 9.40 157.76 194.37 1.00 4511 2829
## age 1.51 0.51 0.49 2.51 1.00 4510 2765
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 47.96 1.76 44.61 51.57 1.00 4301 2784
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).summary(l.ns.hyper)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vk ~ age
## Data: df.ns.large (Number of observations: 400)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 115.91 3.11 109.67 121.96 1.00 3876 2959
## age 4.89 0.10 4.69 5.07 1.00 3978 2932
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 51.21 1.85 47.67 54.95 1.00 3572 3049
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(hypothesis(l.ns.flat, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Flat Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)
plot(hypothesis(l.ns.diffuse, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Diffuse Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)
plot(hypothesis(l.ns.weak, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Weak Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)
plot(hypothesis(l.ns.inf, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Informative Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)
plot(hypothesis(l.ns.hyper, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Hyper Informative Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)
posterior1 <- posterior_samples(l.ns.diffuse, pars = "b_age")[,c(1,2)]
posterior2 <- posterior_samples(l.ns.weak, pars = "b_age")[,c(1,2)]
posterior3 <- posterior_samples(l.ns.inf, pars = "b_age")[,c(1,2)]
posterior1.2.3 <- bind_rows("prior 1" = pivot_longer(posterior1,c(prior_b_age,b_age)),
"prior 2" = pivot_longer(posterior2,c(prior_b_age,b_age)),
"prior 3" = pivot_longer(posterior3,c(prior_b_age,b_age)),
.id = "id")
modelLME <- lm(vk ~ 1 + age, data = df.ns.large)
ggplot(data = posterior1.2.3,
mapping = aes(x = value,
fill = id,
colour = name,
linetype = name,
alpha = name
)) +
geom_density(size = 1.2)+
geom_vline(xintercept = summary(modelLME)$coefficients["age", "Estimate"], #add the frequentist solution too
size = .8, linetype = 2, col = "black")+
scale_x_continuous(limits = c(-20, 20))+
coord_cartesian(ylim = c(0, 1))+
scale_fill_manual(name = "Densities",
values = c("Yellow","darkred","blue" ),
labels = c("diffuse ~ N(5,100) prior",
"weak ~ N(5,10) prior",
"informative ~ N(5, 1) prior") )+
scale_colour_manual(name = 'Posterior/Prior',
values = c("black","red"),
labels = c("posterior", "prior"))+
scale_linetype_manual(name ='Posterior/Prior',
values = c("solid","dotted"),
labels = c("posterior", "prior"))+
scale_alpha_discrete(name = 'Posterior/Prior',
range = c(.7,.3),
labels = c("posterior", "prior"))+
annotate(geom = "text",
x = 0.45, y = -.13,
label = "LME estimate: 0.804",
col = "black",
family = theme_get()$text[["family"]],
size = theme_get()$text[["size"]]/3.5,
fontface="italic")+
labs(title = expression("Influence of Priors on Posterior"),
subtitle = "3 different densities of priors and posteriors and the LME estimate")+
theme_bw()## Warning: Using alpha for a discrete variable is not advised.## Warning: Removed 3653 rows containing non-finite values (stat_density).
and finally again with the large sig data set
summary(mod.lm.s.large)##
## Call:
## lm(formula = vk ~ age, data = df.s.large)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.00463 -0.41362 0.00164 0.45531 1.73849
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.005016 0.032358 0.155 0.877
## age 0.710424 0.034956 20.324 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.647 on 398 degrees of freedom
## Multiple R-squared: 0.5093, Adjusted R-squared: 0.508
## F-statistic: 413 on 1 and 398 DF, p-value: < 2.2e-16summary(l.s.flat)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vk ~ age
## Data: df.s.large (Number of observations: 400)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.01 0.03 -0.06 0.07 1.00 3781 2686
## age 0.71 0.04 0.64 0.78 1.00 3719 2595
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.65 0.02 0.61 0.70 1.00 4063 3084
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).summary(l.s.diffuse)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vk ~ age
## Data: df.s.large (Number of observations: 400)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.00 0.03 -0.06 0.07 1.00 4230 3332
## age 0.71 0.03 0.64 0.78 1.00 4620 2875
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.65 0.02 0.61 0.69 1.00 3717 2671
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).summary(l.s.weak)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vk ~ age
## Data: df.s.large (Number of observations: 400)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.00 0.03 -0.06 0.07 1.00 3898 2346
## age 0.71 0.04 0.64 0.78 1.00 4315 2864
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.65 0.02 0.61 0.69 1.00 4081 3126
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).summary(l.s.inf)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vk ~ age
## Data: df.s.large (Number of observations: 400)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.01 0.03 -0.06 0.07 1.00 4351 3140
## age 0.72 0.04 0.65 0.78 1.00 4213 2665
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.65 0.02 0.61 0.70 1.00 3961 2927
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).summary(l.s.hyper)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: vk ~ age
## Data: df.s.large (Number of observations: 400)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.05 0.15 -0.34 0.23 1.00 2428 2131
## age 3.76 0.13 3.51 4.01 1.00 2054 2299
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 2.90 0.15 2.61 3.21 1.00 2042 2295
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(hypothesis(l.s.flat, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Flat Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)
plot(hypothesis(l.s.diffuse, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Diffuse Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)
plot(hypothesis(l.s.weak, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Weak Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)
plot(hypothesis(l.ns.inf, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Informative Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)
plot(hypothesis(l.s.hyper, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Hyper Informative Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)
posterior1 <- posterior_samples(l.s.diffuse, pars = "b_age")[,c(1,2)]
posterior2 <- posterior_samples(l.s.weak, pars = "b_age")[,c(1,2)]
posterior3 <- posterior_samples(l.s.inf, pars = "b_age")[,c(1,2)]
posterior1.2.3 <- bind_rows("prior 1" = pivot_longer(posterior1,c(prior_b_age,b_age)),
"prior 2" = pivot_longer(posterior2,c(prior_b_age,b_age)),
"prior 3" = pivot_longer(posterior3,c(prior_b_age,b_age)),
.id = "id")
modelLME <- lm(vk ~ 1 + age, data = df.s.large)
ggplot(data = posterior1.2.3,
mapping = aes(x = value,
fill = id,
colour = name,
linetype = name,
alpha = name
)) +
geom_density(size = 1.2)+
geom_vline(xintercept = summary(modelLME)$coefficients["age", "Estimate"], #add the frequentist solution too
size = .8, linetype = 2, col = "black")+
scale_x_continuous(limits = c(-20, 20))+
coord_cartesian(ylim = c(0, 1))+
scale_fill_manual(name = "Densities",
values = c("Yellow","darkred","blue" ),
labels = c("diffuse ~ N(5,100) prior",
"weak ~ N(5,10) prior",
"informative ~ N(5, 1) prior") )+
scale_colour_manual(name = 'Posterior/Prior',
values = c("black","red"),
labels = c("posterior", "prior"))+
scale_linetype_manual(name ='Posterior/Prior',
values = c("solid","dotted"),
labels = c("posterior", "prior"))+
scale_alpha_discrete(name = 'Posterior/Prior',
range = c(.7,.3),
labels = c("posterior", "prior"))+
annotate(geom = "text",
x = 0.45, y = -.13,
label = "LME estimate: 0.804",
col = "black",
family = theme_get()$text[["family"]],
size = theme_get()$text[["size"]]/3.5,
fontface="italic")+
labs(title = expression("Influence of Priors on Posterior"),
subtitle = "3 different densities of priors and posteriors and the LME estimate")+
theme_bw()## Warning: Using alpha for a discrete variable is not advised.## Warning: Removed 3663 rows containing non-finite values (stat_density).
posterior1 <- posterior_samples(ns.s.weak, pars = "b_age")[,c(1,2)]
posterior2 <- posterior_samples(l.ns.weak, pars = "b_age")[,c(1,2)]
posterior3 <- posterior_samples(ns.s.inf, pars = "b_age")[,c(1,2)]
posterior4 <- posterior_samples(l.ns.inf, pars = "b_age")[,c(1,2)]
posterior1.2.3.4 <- bind_rows("prior 1" = pivot_longer(posterior1,c(prior_b_age,b_age)),
"prior 2" = pivot_longer(posterior2,c(prior_b_age,b_age)),
"prior 3" = pivot_longer(posterior3,c(prior_b_age,b_age)),
"prior 4" =pivot_longer(posterior4,c(prior_b_age,b_age)),
.id = "id")
modelLME.l <- lm(vk ~ 1 + age, data = df.ns.large)
modelLME.s <- lm(vk ~ 1 + age, data = df.ns.small)
ggplot(data = posterior1.2.3.4,
mapping = aes(x = value,
fill = id,
colour = name,
linetype = name,
alpha = name
)) +
geom_density(size = 1.2)+
geom_vline(xintercept = summary(modelLME.s)$coefficients["age", "Estimate"], #add the frequentist solution too
size = .8, linetype = 2, col = "grey")+
geom_vline(xintercept = summary(modelLME.l)$coefficients["age", "Estimate"], #add the frequentist solution too
size = .8, linetype = 2, col = "black")+
scale_x_continuous(limits = c(-10, 10))+
coord_cartesian(ylim = c(0, .75))+
scale_fill_manual(name = "Densities",
values = c("blue","red","green","orange" ),
labels = c("small.ns.weak ~ N(5,10) prior",
"large.ns.weak ~ N(5,10) prior",
"small.ns.inf ~ N(5, 1) prior",
"large.ns.inf ~ N(5, 1) prior"))+
scale_colour_manual(name = 'Posterior/Prior',
values = c("black","purple"),
labels = c("posterior", "prior"))+
scale_linetype_manual(name ='Posterior/Prior',
values = c("solid","dotted"),
labels = c("posterior", "prior"))+
scale_alpha_discrete(name = 'Posterior/Prior',
range = c(.7,.3),
labels = c("posterior", "prior"))+
labs(title = "Influence of Weak and Informative Priors",
subtitle = "using a small and large data set")+
theme_bw()## Warning: Using alpha for a discrete variable is not advised.## Warning: Removed 2995 rows containing non-finite values (stat_density).
Below are just null models for each data set. These can be used to get the Bayes factors
Now we’ll compute the actual bayes factors. I’ll walk through two ways to do this but there are others.
The first way involves the library BayesFactor. For this data, we’ll use the function lmBF but you’ll want to explore to see if that is appropriate for your data. lmBF gives us the evidence for the alternative versus a null model. There are other comparisons that can be set up. You’ll need to look into the documentation for your specific needs. A nice feature of lmBF is that it provides an error interval on the BF so you don’t need to run it multiple times to see how reliable it is.
Another way to get a BF is through brms native function bayes_factor. For this, you need to run your null model and then you compare the two. Whichever one goes into the parentheses first is the numerator and the BF you return will be the support for that model. Because the bayes_factor function does not give an error interval, I usually run a few times (here just twice) to ensure the results are equivalent across runs. If they’re not, I increase the iterations in the original model. I then report the average as the BF and the extent of variation as a range or percentage.
lmBF gives the evidence for the alternative. I set bayes_factor to give the evidence for the null. To get values the other way around, just take 1/BF.
#Small n.s. data
bf.s.ns.flat <- lmBF(vk~age,df.ns.small) #.31 (evidence for H1)## Warning: data coerced from tibble to data framebf.s.ns.diffuse <- bayes_factor(null.small.ns.diffuse,ns.s.diffuse) #45.32, 45.93 (evidence for H0)## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5bf.s.ns.weak <- bayes_factor(null.small.ns,ns.s.weak) #5.18, 5.15 (evidence for H0)## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5bf.s.ns.inf <- bayes_factor(null.small.ns,ns.s.inf) #7.23, 7.26 (evidence for H0)## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5bf.s.ns.hyper <- bayes_factor(null.small.ns,ns.s.hyper) #9.56, 9.52 (evidence for H0)## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5#small sig data
bf.s.s.flat <- lmBF(vk~age,df.s.small) #1719.149 (evidence for H1)## Warning: data coerced from tibble to data framebf.s.s.diffuse <- bayes_factor(null.small.s.diffuse,s.s.diffuse) #.04364, .044 (evidence for H0)## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5bf.s.s.weak <- bayes_factor(null.small.s,s.s.weak) #.0048, .00484 (evidence for H0)## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5bf.s.s.inf <- bayes_factor(null.small.s,s.s.inf) #8.67, 8.65 (evidence for H0)## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4bf.s.s.hyper <- bayes_factor(null.small.s,s.s.hyper) #>10000,>10000 (evidence for H0)## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5#large ns data
bf.l.ns.flat <- lmBF(vk~age,df.ns.large) #.12 (evidence for H1)## Warning: data coerced from tibble to data framebf.l.ns.diffuse <- bayes_factor(null.large.ns.diffuse,l.ns.diffuse) #150.58, 150.17 (evidence for H0)## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4bf.l.ns.weak <- bayes_factor(null.large.ns,l.ns.weak) #16.74, 16.72 (evidence for H0)## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5bf.l.ns.inf <- bayes_factor(null.large.ns,l.ns.inf) #6392.59, 6360.84 (evidence for H0)## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5bf.l.ns.hyper <- bayes_factor(null.large.ns,l.ns.hyper) #>10,000, >10,000 (evidence for H0)## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5#Large s data
bf.l.s.flat <- lmBF(vk~age,df.s.large) #5.85 +- .01% (evidence for H1)## Warning: data coerced from tibble to data framebf.l.s.diffuse <- bayes_factor(null.large.s.diffuse,l.s.diffuse) #<.0001,<.0001 (evidence for H0)## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5bf.l.s.weak <- bayes_factor(null.large.s,l.s.weak) #<.0001,<.0001 (evidence for H0)## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5bf.l.s.inf <- bayes_factor(null.large.s,l.s.inf) #<.0001,<.0001 (evidence for H0)## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4bf.l.s.hyper <- bayes_factor(null.large.s,l.s.hyper) #>10000,>10000 (evidence for H0)## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5Now that we’ve looked at the bayes factors for these models if our priors suggest an effect, let’s look at them if our priors suggest there’s not an effect. We’re doing this because the BFs were pretty wonky in some spots above because I set a prior that was ridiculous. When set one that is more reasonable, our results are less wonky. I’ll run all the brms in one block here.
plot(hypothesis(ns.s.diffuseN,"age=0"),plot=FALSE)[[1]]+labs(title="Small, NS, Diffuse Prior")+xlim(-1.5,1.5)
plot(hypothesis(ns.s.weakN,"age=0"),plot=FALSE)[[1]]+labs(title="Small, NS, Weak Prior")+xlim(-1.5,1.5)
plot(hypothesis(ns.s.infN,"age=0"),plot=FALSE)[[1]]+labs(title="Small, NS, Inf Prior")+xlim(-1.5,1.5)
plot(hypothesis(ns.s.hyperN,"age=0"),plot=FALSE)[[1]]+labs(title="Small, NS, Hyper Prior")+xlim(-1.5,1.5)
plot(hypothesis(s.s.diffuseN,"age=0"),plot=FALSE)[[1]]+labs(title="Small, Sig, Diffuse Prior")+xlim(-1.5,1.5)
plot(hypothesis(s.s.weakN,"age=0"),plot=FALSE)[[1]]+labs(title="Small, Sig, Weak Prior")+xlim(-1.5,1.5)
plot(hypothesis(s.s.infN,"age=0"),plot=FALSE)[[1]]+labs(title="Small, Sig, Inf Prior")+xlim(-1.5,1.5)
plot(hypothesis(s.s.hyperN,"age=0"),plot=FALSE)[[1]]+labs(title="Small, Sig, Hyper Prior")+xlim(-1.5,1.5)
plot(hypothesis(l.ns.diffuseN,"age=0"),plot=FALSE)[[1]]+labs(title="Large, NS, Diffuse Prior")+xlim(-1.5,1.5)
plot(hypothesis(l.ns.weakN,"age=0"),plot=FALSE)[[1]]+labs(title="Large, NS, Weak Prior")+xlim(-1.5,1.5)
plot(hypothesis(l.ns.infN,"age=0"),plot=FALSE)[[1]]+labs(title="Large, NS, Inf Prior")+xlim(-1.5,1.5)
plot(hypothesis(l.ns.hyperN,"age=0"),plot=FALSE)[[1]]+labs(title="Large, NS, Hyper Prior")+xlim(-1.5,1.5)
plot(hypothesis(l.s.diffuseN,"age=0"),plot=FALSE)[[1]]+labs(title="Large, Sig, Diffuse Prior")+xlim(-1.5,1.5)
plot(hypothesis(l.s.weakN,"age=0"),plot=FALSE)[[1]]+labs(title="Large, Sig, Weak Prior")+xlim(-1.5,1.5)
plot(hypothesis(l.s.infN,"age=0"),plot=FALSE)[[1]]+labs(title="Large, Sig, Inf Prior")+xlim(-1.5,1.5)
plot(hypothesis(l.s.hyperN,"age=0"),plot=FALSE)[[1]]+labs(title="Large, Sig, Hyper Prior")+xlim(-1.5,1.5)
And finally, compute those BF. the BF on the flat prior won’t change, so I skipped that step.
#Small n.s. data
bf.s.ns.diffuseN <- bayes_factor(null.small.ns.diffuse,ns.s.diffuseN) #45.46, 45.60## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5bf.s.ns.weakN <- bayes_factor(null.small.ns,ns.s.weakN) #5.16, 5.12## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5bf.s.ns.infN <- bayes_factor(null.small.ns,ns.s.infN) #1.09, 1.1 ## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5bf.s.ns.hyperN <- bayes_factor(null.small.ns,ns.s.hyperN) #1,1.001## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5#small sig data
bf.s.s.diffuseN <- bayes_factor(null.small.s.diffuse,s.s.diffuseN) #.04, .044## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4bf.s.s.weakN <- bayes_factor(null.small.s,s.s.weakN) #.0048, .0044## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5bf.s.s.infN <- bayes_factor(null.small.s,s.s.infN) #.0005, .0005## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5bf.s.s.hyperN <- bayes_factor(null.small.s,s.s.hyperN) #.00009, .03## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5#large ns data
bf.l.ns.diffuseN <- bayes_factor(null.large.ns.diffuse,l.ns.diffuseN) #150.48, 150.38## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4bf.l.ns.weakN <- bayes_factor(null.large.ns,l.ns.weakN) #14.98, 15.07## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4bf.l.ns.infN <- bayes_factor(null.large.ns,l.ns.infN) #1.75, 1.79## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4bf.l.ns.hyperN <- bayes_factor(null.large.ns,l.ns.hyperN) #1.01, 1.01## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4#Large s data
bf.l.s.diffuseN <- bayes_factor(null.large.s.diffuse,l.s.diffuseN) #<.0001,<.0001## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5bf.l.s.weakN <- bayes_factor(null.large.s,l.s.weakN) #<.0001,<.0001## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5bf.l.s.infN <- bayes_factor(null.large.s,l.s.infN) #<.0001,<.0001## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5bf.l.s.hyperN <- bayes_factor(null.large.s,l.s.hyperN) #<.0001,<.0001## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5